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ABSTRACT 


A new  method,  based  on  empirical  La  Place  transform, 
was  developed  to  find  asymptotically  efficient  estimates  of 
negative  binomial  distribution  parameters.  These  estimates 
were  found  fairly  close  to  those  found  by  the  method  of 
maximum  likelihood.  Efficiencies  over  95  percent  were 
obtained.  The  method  was  tested  with  a set  of  data, 
generated  by  computer,  and  found  to  be  satisfactory  except 
in  a few  cases.  Maximum  likelihood  also  fails  to  be 
satisfactory  in  these  cases. 
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I.  INTRODUCTION 

The  need  for  readily  computed  parameter  estimates  is 
great.  Maximum  likelihood  estimators  are  known  to  be  asymp- 
totically efficient,  but  in  many  settings  they  are  hard  to 
find.  The  most  popular  alternative  is  the  method  of  moments 
which  usually  yields  readily  computed  estimates,  but  some- 
times these  estimates  are  not  very  efficient.  See  ref.  [2], 

This  study  looks  at  the  efficiency  of  a method  which  uses 
the  probability  generating  function  (empirical  La  Place  trans- 
form evaluated  at  a desirable  value  of  its  argument,  uQ) . 

The  method  presented  herein  requires  computing  power 
greater  than  that  of  the  method  of  moments,  but  less  than 
that  of  maximum  likelihood,  which  calls  for  a psi  function 
capability  . It  can  be  used  with  a hand  held  calculator  that 
has  the  lower  order  transcendental  functions,  i.e.,  logarithm, 
roots  and  power. 

The  basic  idea  is  to  select  a system  of  estimating  equa- 

■ 

tions  which  equates  various  statistics  to  their  expected 
values.  In  chapter  II,  estimating  equations  were  set  and 
efficiency  computations  were  done  using  theoretical  work 
presented  in  [3] . Chapter  III  contains  the  procedure  to  be 
followed, applications  of  the  method  to  some  data  sets,  and 
comparisons  with  the  other  methods. 

In  chapter  IV,  efficiency  and  optimum  uQ  and  t (=-log  uq) 
tables  for  various  p-r  combinations  are  given  and  efficiency 
contours  in  three  different  planes  are  graphed. 
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II.  METHODOLOGY 


J 


A negative  binomial  random  variable  X has  the  probability 
function 


f (X=x;r,p) 


T <r  + x) 
x:r(r)“ 


„x^  r 

q p 


(2.1) 


for  x — 0 , 1 , 2 , ... 
derivatives  of  its 


0<r,  0 < p < 1,  p + q = 1.  The  partial 
logarithm  are,  using  h = log  f 


3h  _ r _ x 

ap  p ~ q 

= ’Mr  + x)  - Y(r)  + log  P 

o IT 


Using  the  basic  recursive  formula  for  the  psi  function  [1] 

x 

4*  (r  + x)  - 'F(r)  = l r J'fVj 

j=l 

one  may  then  express  the  system  of  maximum  likelihood 
equations  as 


x 


0 


Ave  l ^ i + log  P = 0 

1=1, . . . ,n 


(2.3) 
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where  the  values  x^. 


• • • / 


x are  the  data  that  result  when 
n 


the  population  (2.1)  is  sampled  n times.  The  system  (2.3) 
is  nonlinear  in  r,p  and  difficult  to  solve.  Iterative 
methods  must  be  used  and  the  second  member  of  (2.3)  (or 
some  version  thereof)  must  be  recomputed  in  each  cycle. 
This  is  the  main  computational  difficulty  in  using  maximum 
likelihood  in  this  setting. 

The  information  matrix  A,  which  is  defined  as: 


A = -E 


a2h 


32h 


3P 


3p3r 


32h 


32h 


3r3p 


3r 


is  needed  for  efficiency  calculations  and  can  be  developed 
using  the  methods  presented  in  [2] . The  result  is 


A = 


qp 


l 

p 


l 

p 


22 


where 


A22  = r (r)  - E(r  (r  + x)  ) 
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The  determinant  of  A is  found  (using  [2])  as 


|A| 


1 r q11  . r Ini 

2 L (n+1)  ’ (r+n)  I 

^ n=l 


(2.4) 


The  two  estimating  equations  g^  and  g2 > which  are 
exploited  herein,  are  given  as 

g1(r,p;x)  = x - ^ = 0 

(2.5) 

g2(r,p;x)  = tx  - (j^) r = 0 


where 


1 

n 


i=l 


is  the  empirical  generating  function.  (Note:  the  substitution 
t = exp(-u)  converts  this  to  La  Place  transform).  Reduction 
to  a single  equation  is  effected  by  solving  g^  = 0 for  p 
in  terms  of  x and  r. 


p = — (2.6) 

x + r 

and  substituting  this  into  g2-  A function  f (r)  is  obtained 

X — 

in  terms  of  known  variables,  namely  t , x,  and  t. 
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r - tx  = 0 


r + x(l  - t) 


(2.7) 


We  can  solve  equation  (2.7)  for  r by  Newton's  method,  the 
iterative  relationship  being 


f (ri) 

ri+l  = ri  f ' (ri) 


where 


f'(r)  = ( £ )r  • [ (.-  x(i_-  t) } + log( £ )] 

r + x (1  - t)  r + x (1  - t)  r + x (1  - t) 


assuming  an  initial  value  of  r and  a suitable  t can  be 
found  (see  table  4.2).  After  solving  for  r,  p can  be  found 
using  (2.6) . 

Efficiency  for  this  r-p  pair  is  given  by  (see  [2,3]) 


Eff  = 


(2.8) 


where 


jajI 


(2.9) 


| A | is  given  in  (2.4),  the  matices  A and  C in  (2.9)  are 
defined  as 


11 


3p 


12 


r+2 

Eff  = 32 


[£jp)  log(^)  - (a^-)] 2 


3 (l-crt)  ^ r 2 r 2 

r v-tziy r z ^ /i  _+-n 


(1-qt  ) 


r,,  r,.  . 2.  r r q r:n:  , 

-p  (1-qt)  - rqp  (1-t)  ][  I 'TFfnTT1 


(2.12) 


Detailed  computations  are  given  in  Appendix  A.  The  initial 
value  of  r is  found  using  the  method  of  moments  (p,r) , 
and  the  value  of  t used  maximizes  the  efficiency  (2.12) 
when  evaluated  at  p,  r. 


III.  USAGE  AND  APPLICATION 


In  this  chapter  the  procedure  to  be  followed,  in  order 
to  use  the  method  developed,  is  given  in  greater  detail  and 
is  applied  to  three  sets  of  data  taken  from  [4]  and  sixty 
sets  of  data  generated  by  computer  simulation  using  various 
combinations  of  the  parameters. 


A.  PROCEDURE 

The  steps  of  procedure  are  given  below. 

1.  Using  the  method  of  moments,  find  p and  r as 
starting  values.  Starting  values  are  given  by  (3.1) 
(see  [2]) 


~ _ y 

p = x/s 


(3.1) 


r = 


x2 


s2-x 


where  x is  the  sample  mean  and  s is  the  sample  variance. 

2.  Using  p,  r found  in  step  1,  find  the  value  of  t, 
which  maximizes  efficiency  function  (2.12),  and  call  it 


tQ,  and  also  find  the  efficiency,  which  is  maximum  with  this 


pair  p,  r and  tQ,  and  call  it  EFF1 . Tables  (4.1)  and  (4.2) 


have  been  prepared  for  this.  Single  variable  search  and 
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3.  Using  r as  a starting  value  and  tQ,  found  in  step  2, 
find  r*  as  a solution  of  (2.7),  and  p*  using  (2.6).  Newton's 
method  was  used  by  the  author  to  find  r*. 

4.  Step  2 can  be  repeated  with  p*  and  r*  replacing  p 
and  r,  and  using  a new  t (update  t from  Table  (4.2)  if 
the  current  estimated  efficiency  (Table  4.1)  is  not  suffi- 
ciently high) . This  step  involves  the  recomputations  of 
tx  and  is  seldom  needed.  The  new  efficiency  using  p*,  r* 
and  updated  tQ  is  called  EFF2 . 

5.  Stop  if  all  the  following  conditions  are  met. 

(1  - EFF1 ) < e]_  , or  |EFF2  - EFF1 1 <_  e2 


and 


|p*  - p|  < e3, 


and 


|r*  - r|  < e4  . 

Then  p*  and  r*  are  the  estimators  of  true  parameters  p and 
r,  and  the  estimated  efficiency  is  EFF2. 

6.  If  one  or  more  of  the  conditions  in  Step  5 are  not 


r* 


EFF1  = EFF2 


and  go  to  step  3 and  repeat  the  steps  following  until  the 
conditions  in  step  5 are  met.  e^,  e3  an<^  e4  are 

stopping  criteria  which  can  be  chosen  by  the  user . 

B.  EXAMPLES 

Examples  one  through  three  are  the  applications  of  the 
method  to  the  Cricket  score  data  of  Reep,  Polard  and 
Benjamin  [4] . Example  four  is  based  on  data  generated  by 
computer  and  contains  sixty  cases. 

Example  1;  Applying  method  developed  to  the  Cowdrey 
data  [4]  the  following  information  is  obtained. 

Sample  mean  = 1.692 

Sample  variance  = 4.343 

p*  = 0.325 

r*  = 0.816 

Optimum  t = 0.459 

Efficiency  = 0.996 

Number  of  Iterations  = 3 

With  the  same  data  maximum  likelihood  estimators  p,  r, 
and  method  of  moment  estimates  p,  r are  given  in  [2]  as 
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0.329 


0.831 


P = 


p = 0.390 


r = 


r = 1.080 


Example  2:  Example  1 was  repeated  for  Barrington  data 
[4]  and  following  information  is  obtained. 

Sample  mean  = 2.095 

Sample  variance  = 4.939 

p*  = 0.346 

r*  = 1.111 


Optimum  tQ  = 0.538 


Efficiency  = 0.995 

Number  of  Iterations  = 3 

p = 0.345  r = 

p = 0.424  r = 


1.014 

1.543 


Example  3:  Like  the  previous  two  examples,  the  following 
information  is  obtained  by  applying  the  method  to  the 
Graveney  data  [4] : 

Sample  mean  = 1.570 

Sample  variance  = 4.474 

p*  = 0.315 

r*  = 0.722 
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Optimum  t = 0.430 

Efficiency  = 0.996 

Number  of  iterations  = 2 

p = 0.317  r = 0.729 

p = 0.351  r = 0.849 

When  the  first  three  examples  are  examined  carefully 
it  is  noticed  that  the  method  developed  in  this  study  is 
almost  as  efficient  as  maximum  likelihood  method.  Estimators 
are  pretty  close  to  maximum  likelihood  estimators  and 
much  better  than  those  found  by  the  method  of  moments. 

Example  4 : This  example  is  based  on  the  data  generated 
by  computer.  Sixty  cases,  each  of  which  has  a different 
sample  size  and  parameters  p and  r.  Four  sample  sizes, 
which  are  15,  30,  50  and  100;  three  r values,  0.5,  2.5  and 
5;  and  five  p values,  0.05,  0.1,  0.3,  0.5,  0.8  are  used. 

In  each  case  the  data  is  generated  for  given  sample  size 
(n) , p and  r.  Results,  which  contain  case  number  (case  no.), 
sample  size  (n) , p,  r,  method  of  moments  estimates  (p,r), 
optimum  uQ  at  p-r  (Initial  uQ) , efficiency  at  the  end  of 
the  first  iteration  (Initial  eff.),  empirical  La  Place 
transform  estimates  (p*,r*),  optimum  uQ  at  p*-r*  (Final  uQ) , 
efficiency  at  p*-r*-final  uQ  (final  eff)  and  number  of 
iterations  (No.  of  itr.),  are  tabulated  in  the  following 
pages.  Data  generated  for  the  sixty  cases  are  given  in 
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Appendix  B.  The  computer  program,  written  in  FORTRAN  IV, 
is  also  given  in  Appendix  C. 


When  all  the  cases  are  examined,  it's  seen  that  esti- 
mates are  not  given  for  cases  2,  3 and  18.  Case  3 provides 

nothing  to  work  with  as  all  samples  are  zero,  yielding 

— 2 2 — 

X = 0,  s =0.  The  other  two  cases  have  s < x which 

signals  trouble  because  the  method  of  moments  estimates 

r < 0,  p > 1;  values  which  cannot  be  used  to  initialize 

our  iteration  scheme  (3.1)  or  the  maximum  likelihood  scheme 

(2.3)  . 

For  each  of  these  cases  the  function  f (r)  of  (2.7) 
is  a non-negative  decreasing  function  of  r which  tends  to 
+oo  as  r -*■  0 and  is  asymptotically  zero  as  r -*•  °°.  Thus  no 
root  exists.  The  maximum  likelihood  approach  yields  a 
comparable  situation. 
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IV.  RESULT  AND  CONCLUSION 


In  this  chapter  efficiencies  for  various  p-r  pairs  were 
given  in  Table  (4.1).  The  method  developed  was  used  to 
calculate  table  entries.  Three  entries  were  given  for  each 
p-r  pair.  Entry  in  the  middle  is  efficiency  value  for  that 
pair  at  optimum  t or  uQ.  Entries  inside  parenthesis  are 
also  efficiencies  and  found  as  follows: 

Let  u = - log  t 
o o 

If  we  miss  uq  by  ±50  percent,  what  would  the  efficiency  be? 
Efficiency  inside  the  parenthesis  above  the  optimum  efficiency 
answers  the  question  for  u^,  = 0.5  uQ , and  other  one  below 
the  optimum  answers  for  - 1.5  uq. 

The  change  u = -log  t (or  rather  its  inverse  t = e u) 
changes  the  probability  generating  tx  function  into  the 
empirical  La  Place  transform.  It  is  seen  that  a 50  percent 
error  in  its  argument  u results  in  only  a minor  degradation 
of  efficiency,  whereas  a 50  percent  error  in  t would  have 
disasterous  effects.  In  all  of  the  examples,  the  initial 
uQ  was  well  within  this  50  percent  range,  the  worse  case 
being  38  percent. 

After  examining  Table  (4.1)  one  can  conclude  that  effi- 
ciency is  monotonically  increasing  with  increasing  p and 
increasing  r for  r greater  than  one  and  minimum  efficiency 
at  optimum  t or  uQ  is  greater  than  0.95. 


r r i 

Following  Table  (4.1)  optimum  tQ  and  uQ  (in  parenthesis) 
were  tabulated  in  Table  (4.2).  One  may  observe  that  t 

o 

decreases  with  p and  increases  with  r,  whereas  u behaves 

o 

oppositely . 

After  Table  (4.2)  efficiency  contours  were  graphed  in 
p-r  plane  (Figure  4.1),  in  mean-r  plane  (Figure  4.2)  and 
in  mean-exp (r)  plane  (Figure  4.3).  Compared  with  the 
Figure  in  [5]  Figure  (4.3)  always  gives  better  efficiency 
values  than  all  the  other  alternatives. 
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TABLE  4.2  OPTIMUM  t AND  u FOR  p*  - r*  (Cont'd) 


MEAN 

Figure  4.3  Efficiency  contours  in  MEAN-EXP (r)  plane 


APPENDIX  A 


» 


X . 


Since  CovU^t  3 ) = 0 for  i t j because  of  independence, 


it  follows  that 


n 


x . 


C12  ° C21 


- I Cov (x. r t 1) 
n L i 

i=l 


Cov (x, tX)  — E(x*tX)  - E(x)E(t  ) 


One  can  calculate  E(x-tx)  using  the  following  relation. 


x-1, 


G(t)  = E(tX)  G'(t)  = E(xt  ) 


Then 


E(xtX)  = tG'(t) 


where 


G(t) 


= (_JE — )r 

' l~qt ; 


taking  the  derivative  and  multiplying  by  t,  E(xtx)  is  found 


as 


E(xtx)  = 


(1-qt) 


44 


Then  it  yields 


'12 


c . rqp^1!!-!) 

21  (l-qt,r+1 


c22  * " E(92>  * 


Et(tx  - (r%t>r>2) 


n 


x. 


= n Var (t) 


= n Var(i  I t x) 


i=*l 


n 


x . 


= | Var(  l t 1) 

i=l 


= Var ( tX)  = E(t2x)  - [ E ( tX ) ] 2 
= G(t2)  - [G ( t) ] 2 


= (- 


- ( 


2r 


1-qt 


i~qt 


As  the  reader  will  recall,  |A|  was  given  by  (2.4)  and  can  be 


written  as  follows. 


|A| 


m r» 

1 f y q 
~ 2 i n+I 
P 


n=l 


r !n! 
(n+r)  ! 


I 

n=m+l 


n+1 


r In!  i 
(n+r) ! J 


Call  the  summation  from  n = m+1  to  ® the  remainder . The 
remainder  can  be  adjusted  such  that  it  will  be  smaller  than 
any  specified  small  value  (e) . 
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I St 


n=m+l 


n rtn!  r 2 

+T  (n+r) i - L n+1 

n=m+l 


Because  <_  1 (reciprocal  of  binomial  coefficient) 


y 3__.  < 

1 n+1  — 

n=m+l 


m+1  _ 

*=*  l ^ 


A value  of  m can  be  found  such  that 


(m+2)p  - e 


< e 


holds.  Then  this  m guarantees  the  remainder  will  be  smaller 
than  e,  and  the  efficiency  equation  of  (2.12)  becomes 


Eff  = 


[i=3£  - aii^-12 


(1+qtV 


v •>  r 2.  . ™ qn  rln!  x 

- pr(l-qt)2-rqp  (1-t)  ](  I J+T  '(r+n)T} 

n*l 
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APPENDIX  B 

SAMPLES  GENERATED  BY  COMPUTER 

Samples  generated  by  computer  for  sixty  cases  used  in 
chapter  three,  example  four,  are  given  in  the  following 
pages  case  by  case. 

See  subroutine  GNRNB  of  computer  program-3  of  Appendix 
C for  the  simulation  procedure,  which  was  taken  from  [6] . 
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APPENDIX  C 


COMPUTER  PROGRAMS 

Computer  programs,  to  prepare  Table  (4.1)  and  Table  (4.2) 
(Computer  program  1) , to  find  the  estimates  for  the  examples 
one  through  three  of  chatper  three  (computer  program  2) , 
and  to  generate  the  data  for  the  sixty  cases  of  the  example 
four  of  chapter  three  and  to  find  the  estimates  for  them 
(computer  program  3),  are  given  in  the  following  pages. 

There  are  two  subroutines,  namely  GOLDEN  (Golden  section 
search),  SVS  (Single  variable  search),  and  one  function, 
namely  FCN  (to  evaluate  efficiency  function) , which  are 
common  to  the  three  programs  mentioned  above.  These  are 
given  at  the  beginning  under  the  title  COMMON  SUBROUTINES. 
After  that  the  three  programs  are  given. 

One  who  wants  to  use  any  one  of  the  programs,  must  include 
the  common  subroutines  after  the  main  program  of  the  computer 
program  he  chooses. 
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COMMON  SUBROUTINES 


r 


r 


SUBROUTINE  SVS (FCN , TNOT) 

IMPLICIT  REAL *8  (A-H,0-Z) 

COMMON  P , R, DETLAM 
Q=l. ODO-P 

UPPER=1. ODO/DABS (Q) 

ALOW=-l  .ODO/DSQRT  (DABS  (Q)  ) 

RIOU=0 . 000 00 IDO 
START=0. 9999 
U= START 
S=0 . 0001D0 
A=FCN (U) 

U1=U+S 

B=FCN(U1) 

IF  (B.GT.A)  GO  TO  5 
IF  (B.EQ.A)  GO  TO  6 
M=3 

1 U2=START+S*2 . 0D0** (M-2) 

IF  (U2.GE. UPPER)  GO  TO  99 
A=B 

B=FCN(U2) 

IF (A.EQ.B)  GO  TO  7 

IF  (A.LT.B)  GO  TO  8 

M=M+1 

U=U1 

U1=U2 

GO  TO  1 

7 CALL  GOLDEN (FCN, U1,U2, RIOU, 50, TNOT) 

GO  TO  999 

99  CALL  GOLDEN ( FCN, Ul, UPPER, RIOU, 50, TNOT) 

GO  TO  999 

8 CALL  GOLDEN (FCN, U,U2, RIOU, 50, TNOT) 

GO  TO  999 

6 M=2 

CALL  GOLDEN (FCN,U,U1,RIOU,50, TNOT ) 

GO  TO  999 
5 M=3 

2 U2=START-S  * 2 . 0D0  * * (M- 3 ) 

IF  (U2.LE.ALOW)  GO  TO  77 
B=A 

A=FCN(U2) 

IF  (A.EQ.B)  GO  TO  17 

IF  (A.GT.B)  GO  TO  18 

U1=U 

U=U2 

M=M+1 

GO  TO  2 

17  CALL  GOLDEN (FCN, U2,U, RIOU, 50, TNOT) 

GO  TO  999 

18  CALL  GOLDEN (FCN, U2,U1, RIOU, 50, TNOT) 

GO  TO  999 

77  CALL  GOLDEN (FCN, ALOW, U, RIOU, 50, TNOT) 
999  RETURN 

END 
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SUBROUTI NE  GOLDEN (FCN , ENDL , ENDR , RIOU , NMAX , XOLD ) 
IMPLICIT  REAL* 8 (A-H,0-Z) 

COMMON  P , R, DETLAM 

SIGMA= ( 3 . ODO-DSQRT ( 5 . ODO) )/2.0D0 

AL=ENDR-ENDL 

X OLD= END  L+  S I GMA * AL 

FOLD=FCN (XOLD) 

DO  1 1=2, NMAX 

XNE W=ENDL+ENDR- XOLD 

FNEW=FCN (XNEW) 

IF  (FNEW. LE.FOLD)  GO  TO  5 

IF  ( XNEW. LT. XOLD)  GO  TO  6 

ENDR=XNEW 

AL= ENDR- ENDL 

GO  TO  30 

ENDL=XNEW 

GO  TO  20 

IF  ( XNEW. LT. XOLD)  GO  TO  7 

ENDL=XOLD 

XOLD = XNEW 

FOLD=FNEW 

GO  TO  20 

END R= XOLD 

GO  TO  8 

IF (AL.LE.RIOU)  GO  TO  40 

CONTINUE 

RETURN 

END 


FUNCTION  FCN (X) 

IMPLICIT  REAL  *8  (A-H,0-Z) 

COMMON  P,R, DETLAM 
ZR  = 1.0D-73 
Q=1 . 0D0-P 
A=R*P** (R+2.DOJ/Q 
- B=1 . 0D0-Q*X 

C= (DABS ( (B*DLOG (B/P) /P-Q* (1 . 0D0-X) /P ) ) )**2.D0 
D=B** ( 2 . D0*R+2 . ODO) / ( 1 . 0D0-Q*DABS (X) **2.D0)**R 
E=P**R*B**2.D0+R*Q*P**R* (DABS (1.0D0-X) ) **2.D0 
G=-A*C/DETLAM 
F=D-E 

IF  (DABS(F) . LT.ZR)  GO  TO  5 

FCN=G/F 

GO  TO  7 

FCN=-1. ODO 

RETURN 

END 
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COMPUTER  PROGRAM  1 


IMPLICIT  REAL* 8 (A-H,OZ) 

EXTERNAL  FCN 

DIMENSION  PR (14)  ,RR(26)  ,EFF(26,14)  ,TZERO(26,14)  ,TTEST(52)  ,EET’1(26/ 
114) ,EFF2 (26,14) 

COMMON  P,R,DETLAM 
DATA  EFF , TZERO/7 28*0. 0D0/ 

DATA  EFFl,EFF2/728*0. 0D0/ 

DATA  PR/ . 05D0 , . 06D0 , . 07D0 , . 08D0 , . 09D0 , . IDO , . 2D0 , . 3D0 , . 4D0 , . 5D0 , . 6D 
10,.7D0,.8D0,.9D0/ 

DATA  RIV,.5DO,1.DO/1.5DO,2.DO/2.5DO,3.DO,3.5DO,4.DO,4 .5D0,5.D0,6.D0 
1, 7. DO, 8. DO, 9. DO, 10. DO, 11. DO, 12. DO, 13. DO, 14. DO, 1?.  DO, 16. DO, 17 .DO, 18 
l.D0,19.D0,0.0D0,0. 0D0/ 

DATA  TTEST/52*0.0D0/ 

WRITE  (6,444) 

444  FORMAT  ( ' 1 ' ) 

EPS=0.0001D0 
DO  1 J=1 , 14 
P=PR(J) 

Q=1.0DO-P 

K=1 

7 TEST^Q**(K+1)/(P*(DFLQAT(K)+2.0D0)) 

K=K+1 

GO  TO  7 

8 M=K 

WRITE  (6,21)  M 

21  FORMAT (15X,  'UPPER  BOUND='  ,18) 

DO  2 1=1,24 
R=RR(I) 

DETIAM=0 . 0D0 
DO  10  N=1,M 
RA!ITCN=1.0D0 
DO  12  JJ=1,N 

12  RATICN=RATICN*DFLQAT(JJ)/(DFLQAT(JJ)+R) 

DETIA^DETLAM+Q*  *N*RATICN/DFD3AT  (Nfl) 

10  CONTINUE 

CALL  SVS(FCN,TNCT) 

TZERO(I,J)=TNOT 
TPL=DEXP  ( 0 . 5D0*DL0G  (TNQT) ) 

IMUMDEXP  (1 . 5DO*DLOG  (TNOT) ) 

EFF2(I,J)=-FCN(TPL) 

EFF1  (I , J)  =-FCN  (TWIN) 

EFF  (I,  J)  =-FCN  (TOOT) 

2 CONTINUE 

1 CONTINUE 

WRITE  (6,200)  (PR(I) ,1=1,14) 

200  FORMAT  (///,  50X,  'P  VALUES' ,/,12X,14  (2X,F5. 3, IX)  ,/,12X,14  (2X, ' 

1',1X),/) 

DO  13  1=1,24 

WRITE  (6,202)  (EFF1 (I,J) , J=l,14) 

WRITE  (6,201)  RR(I) , (EFF(I,J) ,J=1,14) 

WRITE  (6,202)  (EFF2 (I, J) , J=l,14) 
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201  FOFMAT(3X,  '!$=' ,F5.2,2X,14  (2X,F5.3,1X) ) 
WRITE  (6,203)  (TZERO(I,J) ,J=1,14) 

202  FQFMAT(12X,14 (IX, ' ( ' ,F5.3, 1 ) ' ) ) 

203  FORMftT(5X, 'QP.TNOT' ,14 (2X,F5.3,1X) ,//) 

13  CONTINUE 

STOP 

END 


COMPUTER  PROGRAM  2 


IMPLICIT  REAL* 8 (A-H,0-Z) 

EXTERNAL  FCN 

DIMENSION  FRQ(20) ,OBS(20) 

COMMON  P,R,DETLAM 
WRITE  (6,444) 

444  FCBMAT('l') 

EPS=0.0001D0 
WRITE  (6,91) 

91  FORMAT (5X, 'READ  IN  NO.  OF  OBS.  AND  VARIATES  IN  215  FORMAT’,/) 

READ  (5,81)  NOBS, . NVAR 

81  FORMAT (215) 

WRITE  (6,92) 

92  FORMAT (5X,  'NOW  GIVE  ME  VARIATE  VALUES  AND  CORRESPONDING  FREQUENCY 

UN  THE  FORMAT  2F10.5  FOR  EACH  CARD',/) 

SUM=O.ODO 

SUM1=0 . 0D0 

DO  4 1=1, NVAR 

READ  (5,82)  CBS (I) ,FRQ(I) 

82  FORMAT (2F1 0.5) 

SUM=SUM+OBS (I) *FRQ(I) 

SUM1=SUM1+ (OBS (I) **2 . 0D0) *FRQ (I) 

XBAR=SUM/DFLQAT  (NOBS ) 

SSQR= (SUM1-DFDOAT (NOBS) *XBAR**2 . 0D0) /DFLOAT (NOBS-i) 

WRITE  (6,209)  XBAR,  SSQR 

209  FORMAT (15X,  'USING  GIVEN  DATA  FOLLOWING  INFORMATION  IS  OBTAINED' ,// 

1,20X, 'SAMPLE  MEAN  =' ,F7. 4,//, 2 OX, 'SAMPLE  VARIANCE  =',F7.4,// 

2) 

P=XBAR/SSQR 

R=XBAR**2 . 0D0/  (SSQR- XBAR) 

Q=1 . 0D0-P 
K=1 

7 TEST=Q** (K+l)/ (DFLOAT (K+2) ) 

IF (TEST.LE.EPS)  GO  TO  8 
K=K+1 

GO  TO  7 

8 M=K 
DETLAM=0.0D0 
DO  10  l*=l,M 
RATia*=1.0D0 
DO  12  JJ=1,N 

12  RATION=RATICN*DFLQAT  ( JJ)  / (DFLOAT  ( JJ)  +R) 

DETLAtt=DETTAM+Q*  *N*RATICN/D  FLOAT  (N+l ) 

10  CONTINUE 

L=1 

95  CALL  SVS(FCN,TNOr) 

EFF1=-FCN  (TNCT) 

RSTAR=R 
PSTAR=P 
SUM2=0.0D0 
DO  14  J=1,NVAR 

14  SUM2=SUM2+FRQ  (J)  *TNOT**OBS  ( J) 

TXBAR=SUM2/DFLQAT  (NOBS) 
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CALL  COMPR(TNar,TXBAR,XBAR,RSTAR) 

P=R/  (R*XBAR) 

Q=1.0DO-P 

K=1 

TEST=Q**  ( K+l ) / (P^DFLQAT  (K+2 ) ) 

IF  (TEST. I£. EPS)  GO  TO  88 

K=K+1 

GO  TO  77 

tt=K 

DEHA»=0.0D0 
DO  51  N=1,M 
RATIONAL.  0D0 
DO  52  JJ=1,N 

RATICN=PATION*DFTQAT  ( JJ)  / (DFLOAT  ( JJ)  +R) 
DETLA^DETLAM+Q**N*RATION/DFLOAT  (N+l ) 

EFF2=-PCN(TN0T) 

IF (DABS (P-PSTAR) .GT. 0.005)  GO  TO  29 
IF (DABS (R-RSTAR) .GT. 0.005)  GO  TO  29 
IF  (DABS  (EFF1-EFF2)  .GT. 0.04)  GO  TO  29 
WRITE  (6,94)  P,R,EFF2,TNCT,L 

FORMAT (20X, 'ESTIMATE  OF  P =' ,F7.5,//,20X,  ESTIMATE  OF  R 

17.4, //,20X,EFFICIENCY  =' ,F7. 5,//, 20X, 'OPTIMUM  TNOT 

27.5, //,20X, 'NO. OF  ITERATIONS  =’,FI4,//) 

GO  TO  999 

L=L+1 
GO  TO  95 
STOP 
END 


SUBROUTINE  COMPR(T,TX,XB,RST) 

IMPLICIT  REAL* 8 (A-H,0-Z) 

COMMON  P,R,DETLAM 
R1=RST 

B+R1+XB* (1. ODO-T) 

T¥-  (Rl/B)  **R1 
F=TX-TY 
FP1=-T¥ 

FP2=  (XB*  (1 . ODO-T)  /B)  +DL0G  (Rl/B) 

FP=FP1*FP2 

R2=R1-F/FP 

IF (DABS (R2-R1 ) .LE.0.00001D0)  GO  TO  66 

R1=R2 

GO  TO  30 

R=R2 

RETURN 

END 
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COMPUTER  PROGRAM  3 


IMPLICIT  REAL*8  (A-H,0-Z) 

EXTERNAL  FCN 

DIMENSION  FRQ(200) ,OBS(200) ,PVECT(10) ,KVECT(10) ,NVECT(10) ,DAT(200) 
1,ISD(200) 

COMMON  P,R,DETLAM 

DATA  RVECT/5 . 0D0 , 2 . 5D0 , 0 . 5D0 , 7*0 . 0D0/ 

DATA  PVECT/.8D0, .5D0,.3D0, . IDO , . 05D0 , 5*0 . 0D0/ 

DATA  NVECT/I5 ,30,50,100,6*0/ 

CALL  OJFICM 
WRITE (5, 541) 

541  FORMAT (//,15X,  'READ  ISEED  USED  IN  GENERATING  DATA') 

READ (5,978)  IS 

978  FORMAT (110) 

WRITE (6, 979)  IS 

979  FORMAT (15X, ' ISEED  USED  TO  GENERATE  DATAf' ,110) 

CALL  INr(IS,ISD,200) 

WRITE  (6,444) 

444  FORMAT('l') 

EPS=0. 0001D0 
NCASE=0 
DO  112  IN=1,4 
DO  113  IP=1,5 
DO  114  IR=1, 3 
P=PVECT(IP) 

R=KVECT(IR) 

NOBS=NVECT  (IN) 

NCASE=NCASE+1 
NCP=2*NCASE 
NOHSCP-1 
ISEEDl=ISD  (NCP) 

ISEED2=ISD  (NCM) 

CALL  GNRNB  (NOBS , OBS , P , R,  ISEED2 , ISDEIEDl ) 

CBMX=-99999.D0 
DO  321  LK=l,NOBS 

IF  (CBS(LK)  .GT.CBMX)  OBMX=OBS  (LK) 

321  CONTINUE 

c»iM=cayK+io.DO 

NVAR=0 

390  OBMt4=C»BMX+5.DO 

DO  370  MI=l,NOBS 

IF  (CBS  (Ml)  .IT.OBMN)  OBMN=OBS(MI) 

370  CONTINUE 
ICNT=0 

DO  380  MJ=l,NOBS 
IF  (CBS(MJ)  .NE.OBMN)  GO  TO  380 
ICNT=ICNT+1 
CBS  (MJ)=CBIM 
380  CONTINUE 

NVAR=NVAR+1 

DAT  (NVAR)  =OBMN 

FBQ  (NVAR)  =DFLQAT  (ICNT) 

IF (OBMN.BQ.QBMX)  GO  TO  395 


GO  TO  390 


395  WRITE  (6,889) 

889  FORMAT (//,1  OX, 89  ( '*' ) ,/,40X,  'DATA  GENERATED  IS  GIVEN  BELOW' ,/,39X, 

131 ( '*' ) ,/) 

WRITE  (6,301)  (DAT(LO) ,L0=1,NVAR) 

301  FORMAT(20X, 'VALUES', /,20X, 6 ('-'),/, (5X,20(1X,F5.0) ,/) ) 

WRITE  (6,302)  (FRQ(LO) ,L0=1,NVAR) 

302  fomat920X,  'FREQUENCIES ' ,/, 20X,  11  ( '- ' ) ,/,  (5X,20  (1X,F5. 0) ,/) ) 

WRITE(6,543)  

543  FORMAT (//,10X,  'DATA  WRITTEN  ABOVE  IS  BASED  ON  FOLIOWING  PARAMETERS 
1' ,/,10X,52('-') ) 

WRITE  (6,544)  NCASE,P,R,NOBS,NVAR 

544  FORMAT (/,10X, 'CASE  NO.='  ,15. 5X,  'P='  ,F7.4,5X,  'R='  ,F7.2,5X,  'N=*  ,Ir,5 
IX, 'NO. OF  VALUES=',I3,/) 

SUM=0.0D0 
SUM1=0.0D0 
DO  445  I=1,NVAR 
SUM=SUM+FRQ (I) *DAT (I) 

445  SUM1=SUM1+FRQ  (I) *DAT (I) **2 . 0D0 

XBAR=SUH/DFLQAT  (NCBS) 

SSQR=  (SUM1-DFLQAT  (NOBS)  *XBAR**2. 0D0)  /DFLQAT  (NCBS-1) 

IF  (XBAR.GE.SSQR)  GO  TO  71 
P=XBAR/SSQR 

R=XBAR**2 . 0D0/ (SSQR-XBAR) 

GO  TO  72 

71  R=5.0D0 
P=0 . 9D0 
GO  TO  114 

72  0=1 . ODO-P 
K=1 

7 TEST=Q**  (K+l)  / (DFLQAT  (K+2) ) 

IF  (TEST.LE.EPS)  GO  TO  8 
K=K+1 

GO  TO  7 

8 M=K 
DETLAMHJ.ODO 
DO  10  tt=l,M 
RATICN=1.0D0 
DO  12  JJ=1,N 

12  RATICN=RATION*DET£lAT  ( JJ)  / (DFLOAT  ( JJ)  +R) 

DETLAM=DET1AMK2*  *N*RATION/EFLQAT  (N+l) 

10  CONTINUE 

0=1 

95  CALL  SVS  (FCN,TNOT) 

EFF1=-FCN  (TNCT) 

RSTAR=R 
PSTAR=P 
SUM2=0.0D0 
DO  145  J=1,NVAR 

145  SUM2=SUM2+FPQ(J)  (TNOT**DAT  ( J) 

TXBAR=SUM2/DFIOAT  (NCBS) 

CALL  COMPR(TNOT,TXBAR,XBAR,RSTAR) 
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P=R/  (R+XBAR) 

Q=1.QD0-P 

K=1 

77  TEST=Q**  (K+l)  / (P*DFLQAT  (K+2) ) 

IF  (TEST.IE.EPS)  GO  TO  88 
K=K+1 
GO  TO  77 
88  M=K 

DETLAM=0.0D0 
DO  51  N=l/M 
RATTON=1.0D0 
DO  52  JJ=1,N 

52  RATIC^RATICN*DFIOAT  ( JJ)  / (DFLOAT  ( JJ)  +R) 

51  DETLAtM)ETLAMK)*^RATICN/DFIOAT(N+l) 

EFF2=-FCN  (TNOT) 

IF  (L.NE.l)  GO  TO  141 
UNCT=-DIOG (DABS (TNOT) ) 

WRITE (6, 98)  PSTAR,FSTAR,P,R/EFF1,EFF2,UNOT 
98  FORMAT (/,11X, 57 ('&')  ,/,15X,  'VALUES  OBTAINED  AT  THE  END  OF  FIRST  IT 
1ERATION'  ,/,12X,  'PTML'  ,3X,  'RTHTL'  ,3X,  'PSTAR'  ,3X,  'RSTAR' , 3X,  'ETHTL' 
2,3X,  'ESTAR' ,3X,  'UNOT'  ,/,12X,F5.3,2X,F7.3,2X,F5.3,2X,F7.3,2X,F5.3,3 
3X,F5.3,2X,F7.4,/ , 11X,  57  ( * & * ) ,/) 

141  IF (DABS (P- PSTAR) .CT. 0.001)  GO  TO  29 

IF (DABS (R-RSTAR) .GT. 0.005)  GO  TO  29 
IF (DABS (EFF1-EFF2) .GT.0.01)  GO  TO  29 
UNCT=-DIOG  (DABS  (TNOT) ) 

WRITE  (6,94)  XBAR,SSQR,P,R,EFF2,UNCT,L  

94  FORMAT (15X,  'USING  THE  ABOVE  DATA  FOLLOWING  INFORMATION  IS  OBTAINED 
l/,/,15X,54 ('-' ) ,/,10X,  'SPLE.MEAN'  ,2X,  'SFLE.VAR'  ,3X,  'EST.OF  P'  ,3X, ' 
2EST.OF  R'  ,6X,  'EFF'  ,6X,  'OP.UNOT'  ,3X,  'NO.CF  ITR'  ,/,10X,7  (9  ( '-' ) ,2X) , 
3/,10X,F9.3,2X,F9.3,2X,F9.3,2X,F9.3,2X,F9.3,2X,F9.4,5X,I3,/ /) 

GO  TO  114 
29  L=DKL 

GO  TO  95 
114  CONTINUE 

113  CONTINUE 

112  CONTINUE 

STOP 
END 
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SUBROUTINE  GNRNB(N,OBS,P,R,ISEED2,ISEEDl) 
REAL* 8 P,R,OBS 

DIMENSION  OBS (200) ,GAM(200) ,K(1) 

P1=SNGL (P) 

R1=SNGL(R) 

CALL  GAMA  (R1,ISEED1, GAM, N) 

DO  1 1=1, N 

RIAM=GAM(I)  * (l.-PD/Pl 

CALL  GGPOSH(RLAM,ISEED2,l,K,IER) 

OBS(I)=DBLE(FLOAT(K(l))) 

RETURN 

END 


SUBROUTINE  CQMPR(T,TX,XB,RST) 

IMPLICIT  REAL* 8 (A-H,0-Z) 

COMON  P,R,DETLAM 
R1=RST 

B=R1+XB* (1 . 0D0-T) 

TY=DABS (Rl/B) **Rl 
F=TY-TX 
FP1=TY 

FP2=  (XB*  (1 . ODO-T)  /B)  +DIOG  (DABS  (Rl/B) ) 

FP=FP1*FP2 

R2=R1-F/FP 

IF (DABS (R2-R1) .LE.0.0001D0)  GO  TO  66 

R1=R2 

GO  TO  30 

R=R2 

RETUIN 

END 
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